An Adaptive Finite Element Splitting Method for the 
Incompressible Navier-Stokes Equations 



K. Selim 3 ' 1 ' 2 '*, A. Logg 15 ' 1 - 2 -**, M. G. Larson ' 3 

a Simula Research Laboratory, P.O. Box 134, N-1325 Lysaker, Norway 
b Simula Research Laboratory, P.O. Box 134, N-1325 Lysaker, Norway 
c Department of Mathematics, Umed University, SE-901 87 Umed, Sweden 



Abstract 

We present an adaptive finite element method for the incompressible Navier- 
Stokes equations based on a standard splitting scheme (the incremental pressure 
correction scheme) . The presented method combines the efficiency and simplic- 
ity of a splitting method with the powerful framework offered by the finite cle- 
ment method for error analysis and adaptivity. An a posteriori error estimate 
is derived which expresses the error in a goal functional of interest as a sum 
of contributions from spatial discretization, time discretization and a term that 
measures the deviation of the splitting scheme from a pure Galerkin scheme (the 
computational error) . Numerical examples are presented which demonstrate the 
performance of the adaptive algorithm and high quality efficiency indices. It is 
further demonstrated that the computational error of the Navier-Stokes mo- 
mentum equation is linear in the size of the time step while the computational 
error of the continuity equation is quadratic in the size of the time step. 

Keywords: adaptive finite element method, a posteriori error estimate, 
incompressible Navier-Stokes equations, operator splitting method 



1. Introduction 

Adaptive finite element methods play an increasingly important role in solv- 
ing complex problems in science and engineering. The adaptive methods are in 
general based on a posteriori error estimates, where the error is estimated in 
terms of computable quantities, and adaptive algorithms for automatic tuning 
of critical discretization parameters such as the time step and the local mesh 
size. 
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Several important works have been published on a posteriori error estimates 
for finite element approximations of time-dependent problems, see for exam- 
ple Eriksson et al. [T] , Becker and Rannacher [2] , Giles and Siili [3] , the research 
monograph Estep et al. [J] and the references therein. However, these estimates 
are generally restricted to finite element approximations in space and time and 
thus do not cover commonly used splitting schemes for efficient time stepping. 
Splitting schemes are used to avoid solving coupled systems of equations in 
each time step and have many applications, including reaction-diffusion and 
fluid flow problems. An a posteriori error estimate for a splitting method for 
systems of ordinary differential equations was recently presented by Estep et al. 
[5]. Error analysis of non-Galerkin solutions is considered by Giles and Siili 
[3 with particular focus on error correction; that is, improving the accuracy of 
a computed functional by post-processing. A posteriori error analysis of the 
incompressible Navier-Stokes has been studied in detail before, see for example 
Hoffman 0, Hoffman and Johnson [7], but not for splitting methods. 

In this work, we consider splitting schemes for fluid flow. More precisely, 
we derive an a posteriori error estimate for an incremental pressure correction 
splitting scheme for the incompressible Navier-Stokes equations. This type of 
scheme was originally proposed by Chorin [8] and Temam [9] and was later 
refined by Goda [TO] . Even if a particular scheme is considered, our approach 
extends to other splitting schemes or any other scheme. 

The basic idea is to construct a piecewise polynomial interpolation in time 
of the velocity and pressure and then apply the standard duality argument to 
derive an error representation formula. Since the splitting scheme does not 
satisfy a full Galerkin orthogonality, we are left with an algebraic residual mea- 
suring the effect of the lack of orthogonality caused by the splitting. The final 
estimate thus has three contributions measuring the effect of discretization in 
space, discretization in time and splitting, respectively. A similar approach was 
briefly proposed but not implemented or tested by Bengzon and Larson 
Based on the a posteriori error estimates, we construct an adaptive algorithm 
and investigate the performance of the adaptive algorithm and the quality of 
the error estimate. 

1.1. Outline of this paper 

The outline of this paper is as follows. In the next section, we present 
our model problem (the incompressible Navier-Stokes equations). In Section [3j 
we introduce the inconsistent finite element splitting method that is used to 
solve the Navier-Stokes equations. The a posteriori error analysis is presented 
in Section[4j and the adaptive algorithm is presented in Section[5] The efficiency 
of the adaptive algorithm and the quality of the error estimate is demonstrated 
with numerical examples in Section [6j The paper closes with a summary and 
some concluding remarks in Sectional 
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2. The incompressible Navier Stokes equations 

We consider a fluid governed by the incompressible Navier-Stokes equations. 
For fi C R d (d = 2,3) we seek the velocity u : fi x [0, T] ->• R d and pressure 
p : fi x [0, T] -S- R such that 

u+(u-V)u-V -a(u,p) = / infix (0,T], 

V-it = infix (0,T], ' ^ 

where / is a given body force per unit volume. In ([lj, the first equation is the 
momentum equation and the second equation is the continuity equation. The 
symmetric Cauchy stress tensor a(u,p) is defined as 

a(u,p) = 2ve(u) -pi, (2) 

where v denotes the kinematic viscosity, I is the identity matrix and e(u) is the 
symmetric gradient: 

e(«) = \(Vu+ {Vu) T ). (3) 

We let r = Tq U Tm denote the boundary of fi and associate Dirichlet and 
Neumann boundary conditions with the two disjoint subsets Yd and Tjv, re- 
spectively. On the Dirichlet boundary Yd , we impose a no-slip boundary condi- 
tion (u = 0) and assume a fully developed flow at the Neumann boundary Tjv; 
that is, Vu n — where n is the outward pointing unit normal. This condition 
ensures that the flow does not "creep around the corners" at the inflow and 
outflow. The boundary condition is implemented weakly by dropping the term 
involving V« from the boundary terms, leaving only (f(Vu) T — pi) n. 



3. An inconsistent finite element formulation 

Over the past couple of decades, numerous methods have been developed for 
the numerical solution of the incompressible Navier-Stokes equations. Many of 
these methods are based on a pseudo-compressibility in order to overcome the 
algebraic difficulties of solving the saddle-point problem resulting from a direct 
discretization of the Navier-Stokes equations. A particular type of schemes are 
the so-called splitting schemes where the velocity and pressure variables are 
computed in a sequence of predictor-corrector type steps. 

The first splitting method developed for the Navier-Stokes equations is the 
Chorin projection scheme [HUH]- Chorin's scheme is a so-called non-incremental 
pressure correction scheme where the starting point is to neglect the pressure in 
the momentum equation and solve for a tentative velocity field. The tentative 
velocity is then projected onto a divergence free space, resulting in a Pois- 
son problem for the pressure. In Goda |10j . an improved version of Chorin's 
scheme, the Incremental Pressure Correction Scheme (IPCS), was presented, 
which solves the incompressible Navier-Stokes equations in three steps. In the 
first step, an explicit pressure (the value from the previous time step) is used 
in the momentum equation and in the two subsequent steps, both the pressure 
and the velocity are corrected. 
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In Valen-Sendstad et al. [T2], a comparison is made between different split- 
ting schemes, including a recent scheme by Guermond and Shen [T3], and (sta- 
bilized) Galcrkin finite element methods such as the G2 method by Hoffman 
and Johnson [7] . In this study, six different numerical schemes were were tested 
on six different test problems (making a total of 36 test cases). The test cases 
all involved laminar flow at small to moderate size Reynolds numbers in the 
range 1-1000. For each test problem, convergence in a functional of interest or 
a global error norm was studied for a sequence of refined meshes. The main 
conclusion of this study was that the IPCS scheme was, overall, the most accu- 
rate and efficient method for the particular choice of test problems. Based on 
these results, we choose to base our implementation on the IPCS scheme, but 
emphasize that the analysis is equally valid for any other scheme. 

We consider the IPCS scheme in combination with a Taylor-Hood [14] ap- 
proximation of the velocity and pressure variables; that is, we seek a solution 
U = (U,P) € Vh x Q hl where T4 is the space of continuous piecewise (vector- 
valued) quadratic polynomials and Qh is the space of continuous piecewise linear 
polynomials, respectively. A summary of the IPCS scheme is given in Algo- 
rithm [T] 

To analyze the error of the splitting method, one must construct a suitable 
interpolant/continuous extension in order for the solution to be defined at each 
point (x, t) € f2 x [0, T]. Such an interpolant comes natural for the scheme under 
consideration. Since the solution is computed using a finite element formulation 
in space, we only need to consider interpolation in time. In time, we define the 
discrete solution to be the piecewise linear interpolant on each interval I n based 
on the values U n , U n and P n ~ , P", respectively. For a higher order splitting 
scheme, care must be taken in the construction of the interpolant to maintain 
the order of accuracy. 



4. A posteriori error analysis 



To prove an a posteriori error estimate for the approximate solution of ([I]) 
computed by the inconsistent finite element formulation (the splitting scheme), 
we first state the weak form of ([lj in Section 4.1 and the corresponding weak dual 
problem in Section |4.2| We then derive an error representation in Section |4.3 



from which we obtain the error estimate(s) in Section 4.4 



4-1. The weak primal problem 

The weak form of @ reads: find (u,p) € W = VxQ = {v e L 2 (0,T; [H 1 (Q)] d ) : 
v(-,0) = u ,w| r£ , = 0} x {q e L 2 (n) : (q, 1) = 0} such that 

a((u, P y,(v,q))=L((v,q)) (7) 

for all (v,q) £ W. The test space W is defined analogously to the trial space 
with homogeneous initial conditions. 
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Algorithm 1 The Incremental Pressure Correction Scheme (IPCS) 



Let k n — t n — t n _i denote the time step and /„ = (t n _i,t n ] the corresponding 
time interval. Furthermore, let Vh and Vh denote a pair of trial and test spaces on 
a domain SI. For each time interval I n , we seek the fluid velocity U n = U(-,t n ) G 
Vh and pressure P n = P(-,t n ) G Qh at time t n by solving the following three 
variational problems: 

1) Compute the tentative velocity U* by solving 

<d?(tf*), v) + (a(U n -i,P n - 1 ), e(v)) 
-(v(VU n -i) T n, v) Fn + (P^n, v)r N = (/, v) 

for all v G V/i, including any boundary conditions for the velocity. Here, 
d?(U*) = {U*-U n - 1 )/k n + (U n - 1 -V)U n - 1 andC/™-5 = {U*+U n - 1 )/2. 

2) Compute the corrected pressure P n by solving 

(VP", Vq) = (VP"- 1 , Vg) - fc-^V • U*, q) (5) 
for all q £ Qh, including any boundary conditions for the pressure. 

3) Compute the corrected velocity U n by solving 

(U n , v) = (U*, v) - fc„(V(P" - P"- 1 ), v) (6) 
for all v G Vh, including any boundary conditions for the velocity. 
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The nonlinear form a(-; •) and the linear form L(-) in ^ are defined as 



a((u,p);(v,q)) = / (u + (u ■ V)m, v) + {a[u,p), e(v)) 
Jo 

-(f(Vfi) T n - pn, v) Fn + (V • u, q) dt, (8) 

L((v,q)) = [ T (f,v)dt. (9) 
Jo 

We let r : — > K denote the weak residual of Q ; that is 

r(( M )) = I(( M ))-fl(( Hl p);( M ))= j r\{v,q))dt (10) 



for all (v,<j) € 

4-2. The weak dual problem 

Let now u = (u,p) and let M. = A^(m) denote a given linear goal functional 
(the quantity of interest). The goal functional M. is assumed to be of the form 

M{u) =M T (u(-,T)) + [ M^u^t)) dt. (11) 
Jo 

Here, Ai T and M* describe target functionals at the end time t = T and target 
functionals integrated over the time interval [0,T], respectively. 

We may now introduce the weak dual problem of the incompressible Navier- 
Stokes equations. The dual problem is used below in Section 4.3 to express the 
error in the goal functional A4 in terms of the weak residual (10 1. We let 
z = (z, y) denote the dual solution, where z is the dual velocity and y the dual 
pressure. The (abstract) weak dual problem reads: find z £ W* such that 

a/*(z,v) =M{v) (12) 

for all v £ W*. The dual trial and test spaces are denned by (W*,W*) = 
(W, Wo) where Wo = {v — w : v, w £ W}. The linearized, averaged and adjoint 
form a' : W* x W* — > K is defined by 

a' (v,Su) = a'(5u,v) — / a!(su + (1 — s)U; v)Su ds, (13) 
Jo 

where a 1 denotes the Frechet derivative of the nonlinear form a : W x W —> M 
with respect to its first argument. 

To express the weak dual problem for the incompressible Navier-Stokes equa- 
tions, we start from the abstract dual problem (12 1 and differentiate the non- 
linear form a defined in Q with respect to the velocity field u and the pressure 
field p. The adjoint operator * amounts to replacing the test functions (v, q) 
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in {[sj with the dual variables (z,y), and replacing the linearization variables 
5u = (Su,5p) in (13 1 by the dual test functions (v,q). We find that the dual 
variational problem reads: find (z,y) G W* such that 

/ (z, v) + (z, (u ■ V)v + (v • V)u) + (e(z), a(v, q)) 
Jo 

-(z, u(Vv) T n-qn)r N + (y, V ■ v) dt (14) 
= M T (v(-.T))+ C M\{v,q))dt 



for all (v,q) £ W* , where u = h(U + u). To solve (14), we integrate the first 
term by parts: 

(Z,v)dt=[ (-Z,v)dt + {z(;T),v(;T))-{z(;0),v{;0)). (15) 

Jo 

The boundary term at t = vanishes since (v, q) € W* — W and thus v(-, 0) = 
0. The second term cancels the term Ai T (v(-,T)) in the right-hand side of ( |T4"| 
if we take z(-, T) = ip T where ip T is the (L 2 ) Riesz representer of Ai T . It follows 
that the dual solution may be computed by solving the backward initial value 
problem 

T 

(-£, v) + (z, (u ■ V)v + (v ■ V)u) + (e(z), <r(v, q)) 

-(z, v{Vv) T n-qn) TN + (y, V ■ v) dt (16) 

rT 

M\(v,q))dt, 



with initial condition z(-,T) = i/j t . 

Remark 1. The dual solution may be computed by a direct finite element 
discretization of (16) with u^iU . In particular, it is not necessary to integrate 
by parts the remaining terms of ( 16 ) to move derivatives from the test function v. 

4-. 3. Error representation 

To derive a representation of the error A4(u) — A4(U) = M.{e) in terms of 
the solution z = (z, y) of the dual problem ( [l2| and the weak residual r defined 
in (|lpp , we note that by the definition of the averaged linearized operator a' 
in (13), it follows that 

r r d 

o'(e, v) = / a'(su + (1 — s)U; v)e ds = / — a[su + (1 — s)U; v) ds 

Jo Jo (17) 

= a(u; v) — a(U; v), 
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for all v <E W* , where e = u — U € W* . The error representation now follows 
directly by taking v = e e W* in |l2|): 



M(e) = a' (z, e) = a'(e, z) = a(u; z) — a(U: z) = L(z) — a(u; z) = r(z). (18) 

In other words, the error in the goal function Ai is the (weak) residual of the 
dual solution. 

If now the solution U satisfies the Galerkin orthogonality r(v) — for all 
v S Whk C W, one may subtract a test space interpolant n hk z to obtain M{e) = 
r(z — tt^z) from which the error estimate follows; see Eriksson and Johnson 
[T5] . Becker and Rannacher [5], Bangerth and Rannacher [T5]. However, if the 
solution does not satisfy the Galerkin orthogonality, one must account for the 
lack of orthogonality by adding and subtracting the orthogonality condition. We 
do this in two steps to separately account for the effects of space discretization, 
time discretization and lack of orthogonality: 

■q = M(e) = r(z) 

= r(S - TT h Z + TT h Z - TT hk Z + TT hk z) 

= r(z - ir h S) + r(n h z - ir hk z) + r(-K hk z) 
= Vh + Vk + Vc- 

Here, n h is an interpolant into the semi-discrete space of continuous piecewise 
quadratic vector-valued velocity fields and continuous piecewise linear scalar 
pressure fields at each time t € [0, T], ir k is an interpolant into the semi-discrete 
space of discontinuous piecewise constant functions (at each point x £ f2), and 
TT/ife = T^k^h is an interpolant into the fully discrete test space Whk C W. 

4-4- Error estimates 

To construct an adaptive algorithm based on the error representation (19), 
we estimate rj k and rj c in terms of computable quantities to obtain the total 
error estimate 

\V\ = \Vh + m + T] c \ < \r) h \ + \vk\ + \vd < E h + E k +E C = E, (20) 
where \ij h \ < E h , \r} k \ < E k and \r] c \ < E c . 

The space discretization error estimate E^ 
Starting from the definition rjh = r(z — Tthz), we integrate by parts on each 
cell K g T, where T denotes the triangulation of f2, to obtain 

Vh < [ j2^ Kdt = Eh > ( 21 ) 

J ° KGT 

where r] K = |r^| + {rf^ + \r/ 3 K \ + |^| and 

Vk = (U+(U-V)U-V-a(U,P)-f,z-n h z) K , (22) 

Vk = (M cr ( t/ ' P )1«' z_ Khz)dK\dn, (23) 

Vk = (v^Un, z-ir h z) dKnrN , (24) 

Vk = (V-U,y-7r h y) K . (25) 
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Here, \<j(U,P)\ n = a(U + ,P+)n+ + a(U-,P-)n~ denotes the jump of the 
discrete normal stress a(U,P)n across (interior) edges dK. The time integral 
in (21 1 is evaluated using midpoint quadrature on each time interval /„. 



In practice, we approximate the dual solution z by a numerical approxima- 
tion Z . However, care must be taken when inserting the approximation Z into 
the error representation (19 1 or the error estimate (21 1. In particular, the error 



representation will evaluate to zero if the primal solution satisfies the Galerkin 
orthogonality and the approximate dual solution is computed on the same mesh 
and using the same order as the primal solution. Furthermore, the error esti- 
mate will evaluate to zero since Z — -k^Z = Z — Z = 0. Instead, we compute an 
enhanced version E^Z from the computed dual solution Z by local extrapolation 
on patches, as described in Rognes and Logg [T7]- This allows the dual problem 
to be solved on the same mesh using the same order as the primal problem, 
which has many practical advantages. 

4.4.2. The time discretization error estimate E k 
The time discretization error r) k is estimated by 



\Vk\ = \r(n h z - ir hk z)\ = 
= Ek- 



r*(7TfcS - Tr hk z) dt 



< 



\r l {-K h z - ir hk z)\ dt 
(26) 



To evaluate the estimate E k , we face the problem of integrating the residual 
over each time interval I n . This is challenging as the residual oscillates heavily 
on each interval. For time discretizations defined by a continuous Galerkin finite 
element method in time, the residual is orthogonal to a space of discontinuous 
piecewise polynomial functions on the partition of the time interval [0, T]. As 
a consequence, the residual behaves like a Legendre polynomial on each time 
interval [18) . This is not necessarily the case for a solution computed by a 
splitting method, as is the case here. However, for the sake of analysis, we make 
the assumption that the residual takes its maximum value at the endpoints of 
each interval I n . For a piecewise linear finite element approximation in time, 
the corresponding test space consists of the space of discontinuous piecewise 
constant functions. We may then take the interpolant n k to be the midpoint 
value on each interval to obtain the estimate 



M 



E k < ^fc„|r t (Z(-,t„))-r t ((Z(.,i„_ 1 ) + Z(.,t„))/2)| 



n=l 
M 



(27) 



= \<Tk n \r t (Z{-,t n ))-r t {Z{-,t n ^))\, 



where Z is the approximate numerical solution of the dual problem and M is 
the number of time steps. 



The estimate (27) is used to estimate the size of the time discretization error 



rj k by the adaptive algorithm presented below in Section [5] To control the size 
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of the adaptive time steps, we here derive an alternate estimate. We let i?* 
denote the L 2 Riesz representer of the functional r* and write 

E k = [ \r\n h z - Tx hk l) \ dt = [ \(((n h S - TT hk z, dt 
Jo Jo 

< [ \jn h z - ir hk z 11 Ifli?*! dt 
Jo 

r T (28) 

< max.{kn(t)l^} / K x l% h z - Tt hk z\ dt 

= 5(T)max{fc n (t)|||J?*|||} 
= E k , 

where S{T) = J Q fc" 1 ! 7T/jZ — 7r/jfcz|| dt is a stability factor. The inner product 
(((■, ■))) is here defined by (((u, v))) — (u, v) + (p, q) and the norm ||| • ||| is defined 
as the corresponding norm. We remark that the introduction of inequalities 
in ((28|) may render the estimate less sharp. This however is not a problem 



since (28 ) is not used to estimate the error; it is only used to drive the selection 
of adaptive time steps. 

The norm \\R || of the Riesz representer may be computed directly as follows. 
We first note that the Riesz representer i?* is defined by the variational problem 

v)))=r*{v) (29) 



for all test functions v € Wh- The variational problem (29 1 corresponds to a 
linear system 

M1Z = b (30) 

where M is the mass matrix and 1Z is the vector of degrees of freedom for i?*. 
Clearly, the solution to (30) is given by 1Z = Al~ 1 b. It follows that 

N N N 

i=l j=l »i3=l ' ' 

= {M~ 1 b) T M1Z = b T n, 

by the symmetry of M, The residual norm \R ||| may thus be computed by 
assembling and solving the linear system (30), computing the inner product 
b T TZ and taking the square root. 

4-. 4-3- The computational error estimate E c 

The computational error rj c is computed by a direct evaluation of the weak 
residual for the computed approximate dual solution Z: 

\Vc\ = \r(n hk z)\ w \r(n k Z)\ = \ [ r\n k Z) dt\ = E c , (32) 
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where Z is the approximate solution of the dual problem computed on the 
same mesh and using the same polynomial degree and time steps as the primal 
solution. In our implementation, we have chosen to compute the dual solution 
by a simple application of the dG(0) (backward Euler) method to the linear 
dual problem. We note that, by construction, the computational error estimate 
E c is zero if the primal solution satisfies the Galerkin orthogonality. 



5. Adaptive algorithm 



Based on the a posteriori error estimate derived in Section 4.4 we may now 
formulate an adaptive algorithm for the incompressible Navier-Stokes equations. 
The adaptive algorithm is summarized in Algorithm [2] 

Algorithm 2 Adaptive algorithm 

Given a goal functional M. = M.(u) and a tolerance TOL > 0: 

0) Select an initial coarse mesh and initial time step. 

1) Solve the primal problem ([I]) using (for example) the incremental pres- 
sure correction scheme (Algorithm [lj on the current (fixed) mesh using 
adaptive time steps. 



2) Solve the dual problem (16) backward in time on the same mesh as the 



primal problem and using the same adaptive time steps. 



3) Evaluate the error estimate E — Eh + + E c defined in (21 ), (271 and 



(32), and the error indicators {t]k}- 

4) If E < TOL, then stop. 

5) Refine the mesh based on the error indicators {t)k}- 

6) Continue from step 1). 



In Algorithm [2j we make use of two different tolerances to TOL^ and TOL^ 
which are used to control the errors in the space and time discretization such 
that TOL/j + TOLj. < TOL — E c . The computational error E c is only used as 
part of the stopping criterion E < TOL; it is not used to drive the adaptive 
refinement. However, as will be demonstrated in Section [6] the computational 
error is reduced when the size of the time step is reduced. One may there- 
fore consider extending the adaptive algorithm to control also the size of the 
computational error E c . 

In each adaptive iteration, consisting of a full solution of the primal problem, 
the dual problem and evaluation of the error indicators, the mesh is adaptively 
refined based on fixed fraction marking; that is, a fixed top fraction of the cells 
with the largest indicators are marked for refinement. For mesh refinement, 
we consider two different refinement strategies: the Rivara recursive bisection 
algorithm [TH] and a regular cut algorithm which subdivides all marked triangles 
into four congruent subtriangles and propagates the refinement to neighboring 
triangles using bisection. 
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The step size k n is determined in each time step based on the error estimate 
E k = 5(r)max [ o i T]{^n(i)||| J R*|||}- To achieve E k = TOL fc , we set 

TOL fe TOL fe 

(33) 



S(T)max [tB _ 1)tn] ||J2*|| S(T)lB»V 

where again we have made the assumption that the residual takes its maximum 
value at the endpoints. Since R n is not known until the solution has been 
computed on the time interval /„, which in turn depends on the size of the 



time step k n , it is tempting to replace R n by R 71 in (33). However, this leads 
to oscillations in the time step; if i?™ -1 is large, fc„ will be small and, as a 
consequence, R n will be small, which in turn leads to a large step k n and so on. 
To control the time step, one may introduce a form of smoothing by letting k n 
be the time step determined by 

i. = p^rf (") 

for tolfc = TOL k /S(T) and then take k n to be the harmonic mean 

1k n —\k n 



kn—l 



(35) 



See Soderlind [5D] and Logg [H] for a further discussion on time step selection. 
In practice, we do not compute the stability factor S(T) but instead adjust the 
size of tolfe based on the size of E k . 



6. Numerical results 

We here present numerical results to test the adaptive algorithm and the 
quality of the derived error estimates. An implementation of the adaptive solver, 
including the test problems described in this section, is freely available as part 
of the open source solver package CBC. Solve The package relies on the 

FEniCS/DOLFIN finite element library [23l EH El] . 

6.1. Case 1: Channel flow with wall-mounted body 

As a first test problem, we consider a wall-mounted body (a "flap") immersed 
in a pressure-driven channel flow as illustrated in Figure [l] The kinematic 
viscosity is v = 0.002. As initial condition, we set u = 0. The pressure boundary 
conditions p = 1 at the inflow and p — at the outflow accelerate the flow 
from the initial stationary (zero) solution to the flow field depicted in Figure [2] 
at final time T = 2.5. Note that the solution at final time is not stationary 
which is important when measuring the performance and propagation of time 
discretization errors. 

As a goal functional, we consider the integrated shear stress on the top of 
the flap: 
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u = 



1.0 



p = l 

V« • n = 



1.4 



0.4 



0.6 



p = 
Vm • n = 



u = 



4.0 



Figure 1: Geometry and boundary conditions for the "channel with flap" model problem. 




Figure 2: Fluid velocity (top) and pressure (bottom) at final time T = 2.5 computed with 
fixed time step k = 0.005 and 14 levels of bisection refinement (marking fraction 0.3). The 
final mesh has 16, 581 triangles (76, 085 degrees of freedom). The colorbar indicates the range 
of the scalar pressure field. 
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Figure 3: Dual fluid velocity (top) and dual pressure (bottom) at "final" time t = for the 
channel flow test problem. The colorbar indicates the range of the scalar dual pressure field. 



M\{u) = / / (a(u,p) n) ■ t ds dt 

JO JTi 



\ / o"i 2 (u,p) dsdt = 

JO JFi JO Jrt 



(36) 



viduxfdxi + du2/dxi) ds dt, 



where n = (0,1), t = (1,0) and Ti = {(x 1 ,x 2 ) : 1.4 < xi < 1.8, x 2 = 0.6}. 
As a reference value for the goal functional, we take A4i(u) = 0.0200. This 
reference value was obtained by extrapolation from solutions computed with 
constant time step k = 0.005 on a sequence of adaptively refined meshes. 

6.1.1. Dual solutions and adaptive meshes 

The dual solutions corresponding to the goal functional Ai\ are shown in 
Figure [3] As seen in this Figure, the dual solution clearly reflects the choice of 
goal functional. The dual velocity z (and dual velocity gradients) are large close 
to the top of the flap where the goal functional J\4\ measures the shear stress. 
The same figure displays large spikes in the dual pressure y at the reentrant 
corners. A detail of the dual pressure spikes is displayed in Figure [4] Combined 
with large residuals in the vicinity of the reentrant corners, the result is heavy 
refinement in a region located close to the top of the flap as seen in Figure [5] 
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Figure 4: Detail of the dual pressure field at "final" time t = for the channel flow test 
problem. 




Figure 5: Mesh obtained by 14 levels of recursive bisection refinement with marking fraction 
0.3 for the channel flow test problem. 
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6.1.2. Convergence and efficiency indices 

To investigate the performance of adaptive mesh refinement and the quality 
of computed error estimates, we plot in Figure [6] errors and efficiency indices 
for a sequence of adaptively refined meshes and fixed time step k = 0.005. A 
comparison is made between three different refinement algorithms: recursive 
bisection, regular cut refinement and uniform (non- adaptive) refinement. For 
both recursive bisection and regular cut refinement, we use a fixed fraction 
marking strategy with marking fraction 0.3; that is, in each refinement step, we 
mark for refinement the top 30% of all triangles with the largest error indicators. 

We find that the adaptive algorithm performs very well; a uniformly refined 
mesh requires more than an order of magnitude more degrees of freedom to reach 
a prescribed tolerance. This is evident in Figure [6] by finding the point where 
the error reaches the level |A^i(e)| < 0.001. This level is reached for roughly 
90, 000 degrees of freedom on a uniformly refined mesh, whereas the adaptively 
refined meshes obtained by recursive bisection and regular cut refinement reach 
the same level of accuracy using only 5,000 and 10,000 degrees of freedom, 
respectively. We also note that while the solution obtained by recursive bisec- 
tion converges fastest, the convergence of the solution obtained by regular cut 
refinement is more robust. Computed efficiency indices (error estimate divided 
by actual error) are stable and vary between ca. 3 and 4, which means that we 
overestimate the error, but not by much. 

To study the effect of the choice of marking fraction, we plot in Figure [7] 
errors and efficiency indices for marking fractions 0.1, 0.2, 0.3, 0.4 and 0.5 
for fixed fraction bisection refinement. We note that while a smaller marking 
fraction gives rise to more efficient meshes, that is, a smaller number of degrees 
of freedom are needed to reach a given level of accuracy, more refinement levels 
are needed to reach that level of accuracy. 

6.1.3. Convergence of the global adaptive algorithm 

We next consider the convergence of the global adaptive algorithm. A tol- 
erance TOL = 0.001 is prescribed for the error in the goal functional, here the 



shear stress goal functional M\ defined in (36), and ask the global adaptive al- 
gorithm described in Section [5] to adaptively refine the mesh and select adaptive 
time steps such that |A4i(e)| < TOL. 

Figure [8] shows the convergence of the global adaptive algorithm. The algo- 
rithm converges in five iterations when the prescribed tolerance of TOL = 0.001 
has been reached. Although the actual error reaches the prescribed tolerance 
after only three refinements, the adaptive algorithm performs well; the size of 
the efficiency index is ca. 3. The adaptive time steps are shown in Figure [9] At 
t = 0, the time step is set to the smallest time step from the previous refine- 
ment level. Since the solution is initially at rest, the time residual is initially 
small which leads to an increase in the size of the time steps. As the fluid is 
accelerated by the pressure gradient, the time residual increases and the time 
step is consequently reduced. 



In Figures 10 and 11 we plot the different contributions to the total error 



estimate E = Eh + E^ + E c . Wc find that the error is dominated by the space 
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Figure 6: Error (top) and efficiency indices (bottom) as function of the number of spatial 
degrees of freedom for fixed time step k = 0.005, fixed fraction marking (marking fraction 
0.3) and varying refinement algorithms (recursive bisection, regular cut and uniform) for the 
channel flow test problem. 
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Figure 7: Comparison of errors (top) and efficiency indices (bottom) for varying marking 
fraction using fixed fraction bisection refinement for the channel flow test problem. 
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Error vs error estimate 
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Figure 8: Convergence of the global space— time adaptive algorithm showing errors (top) and 
efficiency indices (bottom) using regular cut refinement with marking fraction 0.3 for the 
channel flow test problem. The given tolerance TOL = 0.001 is reached after five refinements. 
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discretization error Eh, while the time discretization error remains small. This 
indicates that the time steps are unnecessarily small. However, the time steps 
must remain small to preserve stability of the numerical scheme. Although we 
have not taken any special measures to control the size of the time step to 
maintain numerical stability during mesh refinement, the adaptive time step 
selection seems to decrease naturally in each adaptive iteration as a result of an 
increase in the size of the residual |-R§. The computational error E c remains 
practically constant throughout the refinement and we note from Figure [TT] that 
the dominating contribution to the computational error is from the momentum 
equation; the discrete residual of the continuity equation remains small. 



6.1. 4- Convergence as function of h and k 

Finally, we investigate how the error contributions Eh, E^ and E c depend 
on the mesh size h and the time step k. We consider the shear stress goal 
functional M.\ defined in (36) computed on a sequence of uniformly refined 
meshes with mesh sizes h — 0.2, h = 0.1, h = 0.05 and h = 0.025, and fixed 
time steps k = 0.01, k = 0.005, k = 0.0025 and k = 0.00125. 

Figure 12 shows the space discretization error Eh as function of the mesh 
size h. The results indicate that the convergence of the error in the goal func- 
tional is linear with respect to the mesh size. This has not been considered in 
detail but we note that for a P2-P1 Taylor-Hood discretization, we expect the 
convergence of the error in the velocity to be h 3 in the mesh size. However, 
as the goal functional M.\ involves the shear stress, the order of convergence 
is reduced to h 2 . The convergence rate is further decreased by the fact that 
the goal functional considers the shear stress on the boundary and as a result 
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Figure 10: Contributions to the total error E from spatial discretization (E^), time dis- 
cretization (E^) and computational (splitting) error (E c ) for the global space— time adaptive 
algorithm using regular cut refinement with marking fraction 0.3 for the channel flow test 
problem. 
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Figure 11: Contributions to the total computational (splitting) error E c from inexact solu- 
tion of the finite element formulation of the momentum equation (i?™ om ) and the continuity 
equation (E^° n ) using regular cut refinement with marking fraction 0.3 for the channel flow 
test problem. 
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Space discretization error 




Figure 12: Space discretization error for the shear stress goal functional Mi defined in |36| 
as function of mesh size h for varying (fixed) time steps. 



of the singularities at the reentrant corners close to the evaluation of the goal 
functional. We further note from this figure that E^ does not depend on the 
size of the time step with one exception; the error goes up on the finest mesh 
for the largest time step k = 0.01, indicating instability of the numerical scheme 
for large relative time steps. 

In Figure |13[ we plot the time discretization error Ek as function of mesh 
size h and time step k, respectively. We conclude that depends only weakly 
on h and that the convergence of Ek is quadratic in the time step k. 



For the computational error E c displayed in Figures 14 and 15 we similarly 
find a weak dependence on the mesh size h. We further note that the contribu- 
tion from the momentum equation is linear in the size of the time step, whereas 
the contribution from the continuity equation is quadratic. Overall, we thus 
find that the order of convergence is linear in the time step as expected. 

6.2. Case 2: Lid-driven cavity 

As a second test problem, we consider the lid-driven cavity problem on the 
unit square (0, 1) x (0, 1). As boundary conditions, we set u = {xi(l — x\), 0) at 
the top of the cavity (xi = 1) with no-slip boundary conditions on the remaining 
boundary for the velocity. We also fix the pressure p — at the bottom of the 
cavity (x2 — 0). This "unphysical" boundary condition for the pressure gives 
rise to (small) gradients in the pressure field in the vicinity of x-2 = 0. A 
better way to ensure solvability of the pressure update step of Algorithm [I] is to 
require J^p dx = 0. However, we have here chosen to use a Dirichlet boundary 



22 




23 



Computational error (momentum) 
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Figure 14: Momentum computational error E™ oul for the shear stress goal functional Mi 
defined in (|36| as function of mesh size h (top) and time step size k (bottom). 
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Computational error (continuity) 
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Figure 15: Continuity computational error E£ on for the shear stress goal functional Aii 
defined in (|36| as function of mesh size h (top) and time step size k (bottom). 
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Figure 16: Velocity field (left) and pressure (right) at final time T - 
test problem. 



1 for the lid-driven cavity 



condition for the pressure, as this is often used in applications and we wish to 
study its effect on mesh refinement. 

We set the kinematic viscosity to v — 1 and run the simulation over the time 
interval [0, 1]. As a goal functional, we consider a Gaussian- weighted average of 
the ^-component of the velocity field: 



M 2 (u) = 
The weight function (f> is chosen as 



U2(x, t) 4>{x) dx dt. 



(xi,x 2 ) = cexp(-((a;i - xi) 2 + (x 2 - a; 2 ) 2 )/2r 2 ), 



(37) 



(38) 



where {x\,x<z) = (0.75,0.75), r — 0.15 and c « 27.571034 is chosen such that 
J n <j)(xi , x 2 ) dx = 1. As a reference value, we take M2(u) = —0.039389. The 



velocity and pressure fields at final time T = 1 are shown in Figure 16 



6.2.1. Dual solutions and adaptive meshes 

The choice of goal functional generates a source located in the top right 
corner (at X\ — x% — 0.75). The dual solution is advected backwards along the 
primal velocity field and the resulting dual velocity field is shown in Figure |17| 
Notice the large secondary vortex in the top right corner and the small secondary 
vortices in the other three corners. The corresponding adaptive mesh is refined 
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Figure 17: Dual fluid velocity field at "final" time t = for the lid-driven cavity test problem. 




Figure 18: Mesh obtained by 12 levels of regular cut refinement with marking fraction 0.3 for 
the lid-driven cavity test problem (left) and a detailed view of the refined mesh in the top 
right corner (right). 



heavily in the top left and right corners (see Figure 18 1, as a result of pressure 
spikes in these corners, but also at the bottom of the cavity as a result of the 
Dirichlet boundary condition used for the pressure. 
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6.2.2. Error and efficiency indices 

Figure 19 shows the error of the goal functional and the corresponding effi- 
ciency indices for a sequence of adaptively refined meshes, using adaptive time- 
stepping on each refined mesh. Two different adaptive refinement algorithms, 
recursive bisection and regular cut refinement, are compared to uniform re- 
finement. Both adaptive algorithms perform significantly better compared to 
uniform refinement. No significant difference can be noted between the two 
adaptive refinement algorithms, other than that recursive bisection requires ap- 
proximately twice the number of refinement levels to reach the same level of 
accuracy as regular cut refinement. We further note that the efficiency indices 
vary between ca. 1 and 10. Interestingly, the efficiency indices for uniform re- 
finement seem to converge towards 1. 



7. Conclusions 

We have presented an adaptive finite element method for the incompress- 
ible Navier-Stokes equations based on a standard splitting scheme (incremental 
pressure correction). By treating the splitting scheme as an approximation of 
a pure Galcrkin finite clement scheme, one may analyze the error as a sum of 
contributions from space discretization, time discretization and a computational 
error that measures the deviation of the splitting scheme from the pure Galerkin 
scheme. Numerical experiments indicate good performance of the adaptive al- 
gorithm and error estimates that closely match the actual error. The proposed 
method may thus serve as an attractive approach to solving the incompressible 
Navier-Stokes equations, combining the efficiency of a simple splitting method 
with the framework of goal-oriented adaptive finite element methods. 

The presented adaptive algorithm can be further improved by extending the 
adaptive time step selection to control the size of the computational error E c . 
It may also be interesting to consider modified splitting schemes to reduce the 
size of the computational error, in particular the size of the discrete momentum 
residual. 
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